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Introduction 


During this 2 years project, we have made great progress in developing the methods required 
to study bone remodeling as a function of time. The following paragraphs summarizes our 
accomplishments. 

1 CT Acquisition Protocols 

During the two year duration of this project, we have performed a number of experiments to 
determine the best protocol for serial imaging of the calcaneus for quantitative applications. The 
experiments consisted of the following: 

Rob Whalen constructed a phantom that models the calcaneus and its surrounding soft tissue. 
We scanned the phantom using two different energies (80kVp, 120kVp) and 3 different intensities 
(80mA, 200mA, 400mA) in both spiral and conventional modes. For each of the 12 combinations 
of parameters, we obtained 30 scans of the phantom at the same location. We determined the noise 
level for each combination by computing the standard deviation of the voxels within the center of 
the calcaneus model across the corresponding 30 scans. The results of these experiments were: 

(a) We have selected Spiral CT over conventional CT for the following reasons: 

(1) Spiral CT significantly reduces the inter-slice registration error that usually occurs in 
Axial CT and therefore results in more accurate surface and volume reconstruction. 

(2) For thin-slice scanning (as we must do in this application), the patient dose in spiral 
CT is lower compared to conventional for the same longitudinal coverage and 
technique (kV, mA) factors. 

(3) the availability of overlapping slices leads to reduced partial volume errors which in 
term improves the accuracy of the surface detection algorithm, crucial for serial 
registration capability. 

(b) We need to scan at 120 kV and 400 mA in order to achieve noise levels low enough to 
achieve the required precision. 

2 Image Registration 

We have developed, implemented, and tested a surface registration system by modifying and 
extending the methods previously developed for multi-modality applications. As the existing 
system was not optimized for high resolution quantitative computed tomography (QCT), we 
developed new methodologies to improve the registration accuracy and computation speed: 

(a) Surface detection algorithm and representation: 

We have developed a robust and accurate 3-D surface detection algorithm. Our 
implementation of the algorithm results in a set of surface points and corresponding 
normals which are stored in a new efficient data structure. The maximum error of our 3-D 
detection algorithm was 0.2 voxel compared to 0.5 voxel for a conventional detection 
algorithm. The new surface representation and corresponding data structure improved 
registration speed by at least 5 times compared to more traditional triangular mesh- 
generating methods. 



(b) Interpolation function; 

We have developed a methodology to select the optimal interpolation function based on the 
noise structure, voxel size and image content. This significantly reduces the computation 
time and interpolation error incurred in resampling the data volume. 

(c) Gold standard: 

We initially tested our registration method using an external reference system and related 
software as a “gold standard.” During the course of this work, we discovered that our 
new methods had residual errors comparable to errors expected from the gold standard 
registration system. So we developed a new technique for registration using frames that 
resulted in much higher accuracy. This method, which can be used for any stereotaxic 
registration application, allows better quantification of the accuracy of frameless 
registration systems such as ours. As before, it requires the attachment of an external 
reference system that must be kept in place during and between scanning. However, 
existing methods that are widely deployed in neurosurgical applications require accurately 
and precisely constructed frames. In contrast, our method assumes only linear structures 
in the frame and utilizes all 3D information. 

We tested our new “gold standard” using a frame constructed by Robert Whalen at NASA. 
Registration between pairs of successive scans had a maximum error of 0.1 voxel and an 
average error of 0.07 voxel [1][2]. 

We implemented the registration system, including items (a-b) above, together with a user 
interface on a general purpose workstation. This system requires minimal user intervention, 
registers CT volumes containing the whole calcaneus (time required to register two 512 x 512 x 80 
slice volumes is approx. 2 hr., on a SparcStation 20) and is now available for routine use in our 
labs. We scanned 3 calcanei (attached to the Whalen-stereotaxic frame) in different orientations and 
registered each pair (total of 5 pairs) using the new registration system. Based on the newly 
developed gold standard, the registration had an average error of 0.7 voxel (0.18mm) and a 
maximum error of 1.4 voxel (0.36mm). This result clearly showed that the accuracy of the new 
system is suitable for high resolution QCT [3]. 

3 CT attenuation accuracy 

Beam hardening is caused by filtering of the polychromatic x-ray spectrum by the objects in the 
scan field. This effect, if not corrected, can cause severe artifacts in CT images and result in 
inaccurate measures of attenuation. Because we expect bone density in our patient population to 
change over time, beam filtration will change as well. Therefore, the errors caused by beam 
filtration will also change, making it difficult to compare bone density from measurements made at 
different times. We studied the effect that bone density changes have on these errors by scanning a 
specially constructed phantom consisting of a water-filled cylinder containing in its center a thin 
aluminum shell filled with solutions of various concentrations of potassium phosphate to simulate 
different bone densities. We used four different concentrations and 3 different correction 
techniques with the following results: 

(a) Scans using only the standard water correction: Water correction is standard on all 
CT scanners. However, using it alone resulted in errors caused by beam hardening 
that are much greater than our accuracy requirements demand. 



(b) Scans using the iterative bone option QBOl corre ction method: This method is a 
proprietary algorithm available on General Electric CT scanners that corrects for the 
spectral hardening caused by the skull in CT scans of the head. The errors caused by 
beam hardening vary with bone density less than without IBO; however, the errors 
are still too high on the average. 

(c) Scans using a water bolus (with and without IBO); The water bolus pre-hardens the 
beam before it reaches the phantom to reduce beam hardening error. However, the 
corresponding increase in noise level reduced the precision of our measurements 
below our requirements. 


4. Publications 

Reprints or pre-prints of references below are included in the Appendix. 
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Abstract 

This paper presents a new algorithm for frame registration. Our algorithm requires only that the frame be com- 
prised of straight rods, as opposed to the N structures or an accurate frame model required by existing algorithms. 
The algorithm utilizes the full 3D information in the frame as well as a least square weighting scheme to achieve 
highly accurate registration. We use simulated CT data to assess the accuracy of our algorithm. Experimental 
results show that the proposed algorithm is comparable to the best existing technique with knowledge of the 
exact frame model. In situations where there is a discrepancy of more than 0.1cm (0.2% of the frame dimension) 
between the frame and the mathematical model, the proposed technique produces registrations that are 2 times 
more accurate than that from existing techniques. We also outline a routine for estimating the registration error. 
This allows us to take the uncertainty of the frame registration into consideration during an accuracy evaluation of 
a frameless registration system. In one test of using the proposed frame registration system as the gold standard 
to assess the accuracy of a surface-based (frame-less) registration system, error estimates changed by 10% when 
the estimates of the error in the frame-based system were accounted for. 


1 Introduction 

Recent advances in medical imaging have led to an increased interest in image registration. Its usage includes 
stereotaxic surgery, clinical diagnosis, therapy planning, image guided surgery, and many others [1] [2] [3]. Refer 
to [4] [5] for extensive reviews on registration algorithms. 

This work is motivated by our interest in the non-invasive determination of changes in bone mass. Non-invasive 
measurement of changes in bone apparent density require imaging the same skeletal site at different times and 
registering the serial images. To this end, we have developed a semi-automatic, 3D surface-based registration 
technique that does not require the use of an external frame or fiducial markers, we can make non-invasive 
measurements of changes in bone apparent density. Mis-registration can introduce errors and also limit the 
size of the measured volume (the effect of mis-registration errors increases as the size of the measured volume 
decreases). For these reasons, precise evaluation of the registration algorithm accuracy is crucial [6]. 

The availability of a good method for quantifying registration accuracy (Gold Standard) not only allows researchers 
to compare their algorithms, but most importantly, it enables system designers to put registration algorithms 
into practiced use with confidence. Current quantification methods can be divided into three major groups: use of 
1) computer generated data sets; 2) anatomical landmarks and 3) fiducial markers or frame. In the first method 
[7] [8], a second image set is generated from an acquired image set by rotating and translating it through a 
known amount. In this case, the true transformation is known. The drawbacks of this technique are that the 
second image set has (1) noise correlated with the noise in the first set, and (2) lower spatial resolution than 
the first due to the interpolation required to generate it. These effects result in conditions that would not be 



encountered in actual registration of separately acquired scans. The second method utilizes internal landmarks 
to evaluate registration algorithms [9] [10]. Registration error is taken as the root-mean-square (rms) distance 
between corresponding landmarks and is assumed uniform throughout the registered volume. The small number 
of usable landmarks and the errors associated with landmark detection are the main deficiencies of this method. 
The last and most popular method [2] [11] [12] [13] [14] [15] [16] [17], assumes that the transformation of externally 
applied frame/markers is the true transformation. As opposed to anatomical landmarks, the frame/ markers are 
usually constructed in a manner to enable accurate detection. 

The gold standard transformation in categories (2) and (3) by itself has registration error associated with it, 
and that error is commonly ignored. If the error in the gold standard is comparable to that of the algorithm 
under examination, then the evaluation result is useless. This problem is crucial at locations far away from 
landmarks /mar kers/frame as the registration error of the gold standard is likely to change as we move away 
from these structures (usually the landmarks/markers/frame are not in the VOI). In [11] [12], the authors have 
tried to address this problem by inserting rods into the VOI (a cadaver brain) and using the distance between 
corresponding rods as a measure of registration accuracy. This method has good results; however the detection 
error of the rods are ignored and the quantification is precise only in the region near the rods (not the whole 
VOI). Also this method is invasive and can only be performed on cadavers. In general, current quantification 
techniques have two shortcomings: 1) the registration error of the gold standard is ignored, 2) the evaluation is 
accurate only in the region near the landmarks/markers/frame which does not contain the whole VOL 

The current techniques for frame registration are mainly designed for the field of stereotaxic surgery [13] [14] [15] 
[18] [19]. In all these techniques, the authors have assumed that the frame is comprised of parallel N structured 
rods and that an accurate mathematical model of the frame construction is available. Three main shortcomings 
of these techniques are: 1) registration accuracy depends on the accuracy of the frame construction, 2) they do 
not utilize the full 3D information in the frame, and 3) they are not generic; ie. they are not applicable to other 
non-standard frames. 

This paper presents a new algorithm for frame registration. Our algorithm only requires that the frame be 
comprised of a set of straight rods (N structure or exact model are not required), thereby making the algorithm 
applicable to all type of frames. The algorithm utilizes the full 3D information in the frame and a least square 
weighting scheme to achieve highly accurate registration. We also present a routine for estimating the registration 
error. 

Section 2 contains an outline of the proposed registration algorithm and the least square weighting scheme. 
In addition, we define the error measure for registration algorithms comparison and also outline a method for 
estimating the registration error of our proposed algorithm. Section 3 describes the methods we used to validate 
our technique. Section 4 contains the report of the the experimental results from Section 3. Finally, the conclusions 
of our study are presented in Section 5. 

2 Registration Algorithm 

2.1 Problem Definition 

We are given two 3D image sets of a frame. The imaging modality for these images can be the same (acquired 
at different times) or different (e.g. CT and MR). Our goal is to find a transformation (based on the frame) that 
will relate the coordinate systems of these image sets. Before we are able to perform any registration, we need to 
extract the line parameters representing the rods from each image set. There are numerous ways of estimating 
the line parameters; the technique we used is outlined in Section 3.1.3. 

Given the equation of a line L: 


L : r = A;m + a 



where m and a are the gradient and intercept. The equation of a transformed line can be written as 


V : r = fcRm + Ra 4 * t 

where the transformation is comprised of a rotation R followed by a translation t. Note that any rigid body 
transformation can be expressed in terms of R and t. Next, we represent a set of l lines by (M, A, fi m , ft a ) by 


M = 

A = (ai,a 2 , * * * ,aj) 

” (5j m i , 2 m2 , • * • , Smi) 

= (S a i , E a2 , • • • , E a /) 

where m* and a* are the gradient and intercept of the ith line. The estimation error of m< and a* are captured 
in E m i and E ai respectively with 


Emi = £{(m< - E(mi))(ui{ - E( mi)) T } 
E oi = £{(a i - J B(a i ))(a j -E(a i )) T } 


We further assume that the estimation errors for each line are unbiased and independent of those for all other 
lines. So mathematically, the problem can be restated as: 

• given two sets of lines (Mi, and (M 2 , A 2 ,fi m> 2,^a,2) 

• find a transformation (R*,t*) that optimally matches lines in the first set to their corresponding lines in 
the second set. 

2.2 Determining the Rotation and Translation 

We divide the transformation determination problem into two subproblems: 

• Determine the rotation R* using the following criterion: 

i 

R* = argmin || mi ( i — Rmj (2 H2 subject to R r R = I ( 1 ) 

R i=i 

where m ij is the ith column of Mj. 

• Determine then the translation under the following scenario: 

We perform a transformation of (R*, t) on (M 2 , A 2 , ft a(2 ) to obtain (M^, A' 2 (t), For 

each line L* in (Mi, Ai,ft m ,i,ft a ,i) and its corresponding pair LJ(t) from (M' 2 , A' 2 (t ),n' mt 2,n , fl>2 ), we 
define a plane 7 r< which passes through the centroid of the VOI and is perpendicular to L iy as shown in 
Figure 1 . We define p* as the intersection point of and the plane tt<. Similarly, p '*(*) is the intersection 
point of LJ(t) and 7r*. The optimal translation t* is determined by the following criterion: 

t* = argmin]T || p< -p'<(t) ||| . 

*=1 


( 2 ) 




Li 

Figure 1: The optimizing criterion for translation 


The optimizing rotation for Equation (1) is 

R* = UV T (3) 

where U and V are the unitary matrices in the singular value decomposition of MiM^ which is given by 

U t MiM^V = diag(\ u \ 2 ,\3 ) = A 

Refer to Appendix A for a detailed derivation. In Appendix B, we show that Equation (2) is equivalent to a 
quadratic optimization problem and can be solved by using the calculus of variation. The optimizing translation 
is 


= { y ' (I + 7 ~ r, n<nf )} 1 y^( k » m < n . _ 

|=1 ' *— 1 


(4) 


2.3 Weighting Scheme for R* and t* 

In Section 2.2, we have formulated the optimization criteria based on the assumption that the line estimation 
error variances are equal. In cases where f l my and H a are available and the error variances are not constant, we 
could improve the accuracy of the frame registration by using a weighting scheme. In this method, we weight the 
registration error of each line by a weight chosen inversely proportional to their estimation error variances. We 
reformulate Equation (1) and (2) as 

i 

R* = argmin^tu/ii || m^i - Rm^ ill subject to R r R = I (5) 


and 


i 

t* = argmm^tut. || Pi -p'j(t) III 

1=1 


( 6 ) 



with 


E(|| m M - Kmy H*) 

^(ii p» -p', mill) 

where (1R, T) is the true transformation. Appendix C shows the procedure for calculating and xvu . Equation 
(5) and (6) can be solved using the procedures in Appendix A and B. 

2.4 Error Measure 

Before we could perform any comparison between registration systems, we need to define the error measure that 
is to be used in the comparison. Let the true transformation be (#,T). For any other transformation (R,t), the 
registration error at a point p is 


1 

W Ri 

J_ 

Wti 


<* = (*- R)p + (T-t). 

By using the triangle inequality, we have 

E(\\S\\ 2 ) < £(||Rp|| 2 ) + S(||t|| a ) 

< E(\\ max \) HpIIs +®(||t|| 2 ) 

where R = IR-R, t = T- t and |A m **| is the largest absolute eigenvalue of R. E^Xmax]) and £?(|| t || 2 ) are 
good indicators for the rotational and translational registration accuracy respectively. Thus we use (E( 1X^*1) , 
£?(|| t || 2 )) as an error measure for comparing registration algorithms. We denote £(|A mox |) by cr and E(\\ t H 2 ) 
by et . Note that the translational error in a registration is constant throughout the VOI and its expected value 
is equal to e*. Whereas, the rotational error increases linearly with the distance from the center of rotation and 
a multiplicative constant of cr. Thus the overall accuracy of a registration algorithm depends on the value of 
(e*, e t ) and the size of the VOL We refer ( cr , c t ) as the frame registration error with rotational registration error 
cr and translational registration error e*. 

In actual situations, we do not have the luxury of performing many runs of registration to obtain an estimate 
for frame registration error. However, it will be extremely useful for us to know the degree of uncertainty in the 
frame registration. Due to the nice properties of ^(RTR) and E(\\ t H 2 ), we outline their derivations in Appendix 
D. We then approximate cr with the largest eigenvalue of E( R T R)4 and Ct with E(\\ t H!)^- 


3 Experimental Design 

3.1 Registration Accuracy Comparison 

In this section, we compare our technique with two other commonly used techniques. We denote the first technique 
as the Parallel-N technique [13] and the second as the Point-Based technique [15] [18] [19]. The Parallel-N 
technique uses 2 detected N-structures as the coordinate system and the Point-Based technique matches each 2D 
slice separately to an N-bar model. 

3.1.1 Simulation 

Simulated CT data of a frame, shown in Figure 2(a), were obtained by a simulation program [20]. All CT images 
(512x512 pixels) were simulated with the following parameters: Axial, 80kV, 100mA, 0.3x0.3xl.0 cm 3 voxels. 
Figure 2(b) shows a slice from the simulated CT volume. The variance of the poisson noise was experimentally 
adjusted to be similar to that in an actual scan acquired and reconstructed using the identical parameters. 




(a) 


(b) 


Figure 2: (a) The frame model used in the simulation, (b) A slice of the CT scans. 


3.1.2 Data 

Sixty contiguous 1.0cm slices were generated from z — —5cm to z = 55cm with the frame in the orientation 
shown in Figure 2(a). We denoted these sixty slices as a CT volume of the frame. We simulated 50 CT volumes 
separately such that the noise in each CT volume was independent of those in other CT volumes. We then 
transformed the model of the frame by a rotation R=(5° about x, 5° about y and 15° about z) and a translation 
t=(-1.5cm, 1.5cm, -1.5cm). Another set of 50 CT volumes were obtained with the frame in the new orientation. 

3.1.3 Methods 

For every slice, we used an edge detection algorithm to detect the circumference of the rod [12]. We then performed 
a least square fit of the circumference to an ellipse. The center of the ellipse was taken to be the center of the rod. 
Using this method, we estimated the center of the rod in every slice. The parameters of the line representing the 
rod, (m, a, E m ,£ a ) were estimated using the procedures described in Appendices E and F. For each registration 
technique, we performed 50 registrations and compared the results obtained with the known true transformation. 
The exact frame model was used in both the Point-Based and Parallel-N techniques. 

For every registration, we used the procedures in Appendix D to obtain l£(R r R) and E(\\ t ||). We then estimated 
the frame registration error (ci?,ct) and compared it with the true value obtained from the previous experiment. 

To study the sensitivity of the Parallel-N and the proposed technique to model inaccuracies (the frame and its 
model do not match), we repeated the experiment with an inaccurate frame. We formed an inaccurate frame by 
perturbing the end points of all the rods in the frame with a gaussian random variable (with variance a n ). A 
50 CT volume set was obtained using this frame in the orientation shown in Figure 2(a). A the second set was 
obtained with the inaccurate frame in the second orientation, as in Section 3.1.2. We performed 50 registrations 
using each registration technique and compared the results with the true transformation. We performed 50 monte 
carlo simulations for each cr n ranging from (0.0001cm to 0.2cm). 



Method 

Cfl 

e t (cm) 

Point-Based 

Parallel-N 

Proposed 

0.0074 ± 0.000024 
0.0006 ± 0.000004 
0.0009 ± 0.000006 

0.7167 ± 0.00037 
0.0074 ± 0.00007 
0.0053 ± 0.00004 


Table 1: This table shows the expected registration errors (rotation |A max | and translation || t ||) in the various 
techniques, (voxel size: 0.3x0.3xl.0 cm 3 ). 


3.2 Evaluation of the Surface-Based Registration Algorithm 

We used the proposed frame registration system to evaluate the accuracy of the surface-based registration system 
that we have customized from [11] for quantifying bone loss over time. We secured excised human calcanei within 
a custom made fiducial frame shown in Figure 3. We scanned the frarae/calcaneus in different orientations. 
Scanning was done at a resolution of 0.25x0.25x0.5mm 3 . The rods in the image volume were extracted by 
the procedure described in Section 3.1.3. We determined the accuracy of the surface registration algorithm 
by 1) obtaining the frame transformation using our proposed algorithm, 2) obtaining a transformation for the 
calcaneus using surface-based algorithm, 3) transforming all calcaneus voxels using both the frame s and the 
surface’s transformation, 4) calculating the distance (error) between corresponding voxels in the two transformed 
calcaneus images. Note that none of the current frame registration techniques is applicable to this custom-made 
frame. To examine the effects of noise and contrast, we repeated the whole experiment with the frame/calcaneus 
immersed in a water bath. 



Figure 3: An excised human calcaneus tightly secured in the custom made frame. 



4 Results 


4.1 Registration Accuracy Comparison 

For each registration, we obtain (R,t) which is used for calculating (||A ma x||> II t ||). The mean of (||A ma x||> II t ||) 
from the 50 registration runs is used as an estimate for (ch, tt)- The results are tabulated in Table 1. As the 
Point-Based method only utilizes points information on a single slice to perform the registration, its performance 
is not as good when compared to methods that utilize the entire line. The results show that the performance 
of the proposed method is comparable to that of Parallel-N method. The proposed method produces a better 
translational registration while the Parallel-N method is more accurate in rotational registration. The amount 
of information contained in the true frame model and six rods (four vertical, two slanted) enables the Parallel-N 
method to produce a better rotational registration than the proposed method whose registration is solely based 
on the information from eight rods (four vertical, four slanted). However for translational registration, the z-axis 
translation (refer to Figure 2(a)) can only be determined by slanted rods. In the proposed method, the use of 
two additional slanted rods provides substantial amount of information on the z-axis translation which leads to a 
more accurate translational registration in the z-axis. The significant gain in the z-axis translational registration 
accuracy results in a better overall translational registration for the proposed method. For the simulation, the 
mean absolute translational registration error in (x-axis, y-axis, z-axis) for the proposed and the Parallel-N 
methods are (0.0026,0.0021,0.0032) and (0.0020,0.0026,0.0060) respectively. Note that the proposed method is 
two times more accurate than the Parallel-N method in the z-axis translational registration. 

The estimated values of e R and e t are 0.00094 ± 0.000037 and 0.0061 ± 0.00013 respectively. A comparison of 
these figures with Table 1 shows that there is bias of 0.00004 in the estimate of e R and 0.0008 in the estimate 
of e t . Due to the fact that E(X 2 ) > E( \X\) 2 for any random variable X , the approximation done in Section 2.4 
introduces a positive bias into the estimates. As we are going to use these estimates to form an upper bound for 
the registration error of the proposed algorithm, the biases are not significance since they are small (less than 
15% of the true value) and positive. 

Figure 4 shows the performance of the Parallel-N and the proposed techniques on an inaccurate frame. The 
accuracy of the Parallel-N technique is degraded with the increase in frame-model discrepancies. Whereas the 
performance of the proposed technique, which is not based on any frame model, is not affected by the frame-model 
discrepancies. If the mis-match between the frame and the model is more than 0.01cm, the registration error for 
the Parallel-N technique increases exponentially. The performance of the proposed technique surpass that of the 
Parallel-N technique, in both rotational and translational aspect, when the frame-model discrepancies are more 
than 0.05cm (0.1% of the frame dimension). This amount of discrepancies are common in practice which may 
due to the frame being dropped or damaged or mal-machined. 


4.2 Evaluation of the Surface-Based Registration Algorithm 

Figure 5(a) shows an error plot for a slice of the calcaneus (in air) using the frame-based transformation as truth. 
The availability of the estimated (R,t) allows us to account for the frame registration error in our evaluation. 
By the triangle inequality, we have 

registration error \t rue < registration error \ frame 4" frame registration error. 

Thus by adding the estimated frame registration error to the error in Figure 5(a), we have formed an upper 
bound for the true registration error of the surface-based algorithm. Figure 5(b) shows the registration error 
(with respect to the true transformation) of the surface algorithm for that particular slice. 

By repeating the same analysis for the rest of the calcaneus, we have obtained an worst case error map for 
the surface-based technique on the whole calcaneus. Table 2 shows a performance summary for our surface- 
based algorithm. Note that there is a significant increase in registration error (10%) when the frame registration 
uncertainty is taken into consideration. 




(a) (b) 


Figure 4: This figure shows the effects of frame model discrepancies on the performance of the Parallel-N and the 
proposed method, (a) The translation error E(||t||) plotted against the perturbation a n (the standard deviation of 
the gaussian perturbation), (b) The rotation error E(|A ma z|) plotted against the perturbation <r n . The horizontal 
lines are 1 standard deviation from the mean, (voxel size: 0.3x0.3xl.0 cm 3 ) 




Error (mm) 
relative to FYame 


Error (mm) 
true 


Ave I Max I Min I Ave Max 


Calcaneus in Air 
Calcaneus in Water 



Table 2: This table shows the evaluation results for the surface-based algorithm. The first set of results is obtained 
under the assumption that the frame registration represents the true transformation. The second set of results 
accounts for the uncertainty in the frame registration 














(a) (b) 


Figure 5: This figure shows the contour plot of the registration error on a slice (through the centroid of the 
VOI) from the image volume, (a) shows the registration error with the frame’s registration assumed to be the 
true registration, (b) shows the registration error when the uncertainty of the frame registration is taken into 
consideration 

5 Conclusion 

We have described an accurate frame registration algorithm. The proposed algorithm assumes only linear struc- 
tures (rods) in the frame and also utilizes all the 3D information to perform the registration. The proposed 
method is not restricted to the head frame (which is comprised of N-structures). So to perform registration sys- 
tem accuracy evaluations on a VOI other than the head, one only needs to designate some linear structures as the 
frame. A comparison with various current frame registration algorithms shows that the proposed method is more 
robust and more accurate. In addition, we have derived an estimator for the frame registration error. We have 
shown that it is necessary to incorporate frame registration errors into the evaluation of frame-less registration 
system and that doing so gives an upper bound on the expected error of the frame-less system. 
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A Determining the Rotation R* 

Given two matrices Mi and M 2 , we would like to find a rotation matrix R* such that 


R* = argmin V || m<,i - Rm<,2 \\\ subject to R r R = I 

R ■ 1 ^ 

»=1 

The above optimizing criterion can be written as 

R* = arg min || Mj - RM 2 |& subject to R T R = I (7) 

R 

where ||.||f is the Frobenius norm [21]. Since R is orthogonal, then 

|| Mi — RM 2 || trace (Mf Mi) -f trace(M^M 2 ) — 2trace(R T Mi M^) 

then (7) is equivalent to the problem of maximizing trace (R t MiM^). The maximizing R can be found by 
calculating the singular value decomposition (SVD) of MiM^. If 

U t MiM^V = dia^(Ai, A 2 , A3) = A (8) 

is the SVD decomposition, then 

trace(R T Mi M%) = trace(R T UAV T ) 

= trace(V T R T UA) 

= zuAi + Z22 A 2 4* Z33A3 

where Z = V T R T U and is the (i,j) element of Z. As U and V are unitary, the upper-bound is obtained by 
setting Z = I. Thus R* = UV T [21]. 

B Determining the Translation t* 

Let the planes assigned to each line be (7ri , 7r 2 , * • * , 7r/). The point of intersection of a line L : r = a + A;m and a 
plane tt : r.n = d is given by 


d — a.n 

p — m + a. 

m.n 


( 9 ) 


Similarly, the intersection of the plane and a line with a translation t, L:r = a-ft + km, is given by 


, x d - a.n t.n 

p(t) = m -ha m -h t. 

m.n m.n 


So the minimization problem 


t* = argmin || p< - p' f (t) ||^ 


*=1 


can be written as 


t* = arg min V || k< ^-m'j +t ||| 

t ^ m'j.nj 


( 10 ) 



with 


R a i, 2 , 
R*Uli t 2 
di - a 'j.Uj 



m'j.n i 


-m^ 4- a! { - 


at.n i 


mi.ni 


•mi — a* 


where m^ and a* t 2 are the gradient and intercept of the ith line in (M2, A2,ftm,2, ^0,2)- Equation ( 10 ) is a 
quadratic minimization problem and t* can be found by the calculus of variation. At the stationary point, we 

have 




(ni.m'i ) 2 


T 

n i n i 


- — 4— n'jin'fjt* = ‘ ,‘;n< - kj). 

n.m'i “ nj.m * 


(ii) 


With this equation, t* can be found by a matrix inversion. 


C Determining the weights 

Since the errors in m are unbiased and independent, we can show that at the true rotation 

E( m it i - ^mi )2 ) = 0 

£/(|| m< t i — 3£uij f 2 1 1 2 ) = trace (Emi, 1 ) 4* trcice(Ernt,2) 

1 

W Ri 

Similarly, with the assumption that the errors incurred in the line fitting procedure are unbiased and independent, 
at true transformation ($, T) we have 

£(p<-p'<m) = 0 

EdlPi-p'iWWl) = £(||p i -p7ll^) + S(||p' j (0')-p' j (T)||l) 

From ( 9 ), we have 


£(ii p - £(p) iii) = 


SOI a \\l) 

m.n 

*(ll (I - ^)a HI 

m.n 

„._ T/ , mn T ., I nm T 

E a T I )(I )a) 

m.n m.n 


= E( a r 0a) 

= S( 0 )«(S„ M )y 

U 


where© = (I- — )(I - — ) 
v m.n m.n 


Using a similar procedure, we obtain 


ij 

* £(R* T eR-) 0 (£ a< , 2 )y 


E( || p',(f) - p»,(f) ||5) 



Thus, combining these results, we have 

— « V(0)y(SaM) a + £(R ,T eR-)y(5W)y. 

Wti i: a. 

D Derivation of the Registration Error Estimators 

D.l Uncertainty in Rotation 

Without loss in generality, we assume that M 2 is measured exactly and its detection errors are transfered to the 
respective rods in Mj. Based on this assumption, we have 


Mi = SM 2 + Mi - M 2 
= 3tM 2 +A 

where 3? is the true rotation. Our aim is to find a formula for 


E(R t R) = E{( R* - 30 T (R'* - 31)}. (12) 

where R* = UV T . Since we have no access to 3t, we could treat 31 as a small perturbation from R*. Following 
the derivation of R*, let 3? = UV T . It can be shown that U, V,U, V are the eigenvectors of A, B, A + E and 
B-f F respectively, where 

A = MiM^M 2 Mf, B = 

E = AM[M 2 A r - MiM^M 2 A r - AM[M2M[ 

F = M 2 A T AM|’-M 2 MfAMj-M 2 A r MiM?’ 

We can say that U contains the perturbed eigenvectors of A with perturbation E, and V contains the perturbed 
eigenvectors of B with perturbation F. 

The first order eigenvector perturbation [22] of U and V are 


Ujt = u k + 


3 


E 

i=i, i?k 


ufEut 
(A* - X<) 


u», 


"Vfc = v fc + 


3 


E 


ufF v t 
(At — Aj) V< 


(13) 


where the As are the eigenvalues of A. Rewriting (12), 

(R* - 3i) T (R* - 3i) = 21 - VU T 3t - X T UV T 

= 21 - (V + <SV)(U + <SU) T R* - R* t (U + <HJ)(V + <5 V) t 
= -<5V<5U t R* - R* T <SU<fV T - ViU T R* - 5VU T R* - R* t U<TV t - R* T <5UV T 


Thus 

E{{ R* -3?) r (R* -3?)) = -E(6V6U T )R* -R mT E(6XJ6V T ) 

- VE((5U t )R* - E(<TV)U t R* - R* t U£(£V t ) - R* T £(<5U)V T 

where J9(<5V<SU T ), E(SJJ) and E(6V) can be found by substituting E{ E) and E{ F) into (13). 



D.2 Uncertainty in Translation 

FYom Equation (6) and (11) we have 


y> ti (i + - =25- - a=L)i = £t*.(i - - irS")( a < - *'<> 

“ (ni.mt ) 2 Hi. mi n.mj " ni.mi n^nii 


T 

min/ 




mm; 


min; 


which can be written as 


i 

t = $ -1 £ ®<( 5 * - a '»)> 

<=i 


where 


Hence, 

E(t T t) = «C((i4-^) r or)*- r *- 1 Deift -i'O)] 

t=l *=1 

= E £[(a< - a '<) T ©f $ _T $~ 1 0i(a 4 - a',-)]. 

t~l 

Since the errors are independent, we have 


E(i T t) = 53{53(ef*- T $- 1 e,) M (E. il i) M + £(R* T ef*- T *- 1 e i R*W5WW 

«=l pg pg 

E Line Fitting 

A line can be represented by the following equation 


$ 

0j = w t i(I - 


5> ( i +7 i^ 

££ (ni.mi) 2 


T 

min/ 


ni.mi n.mj 


?-) 


Him; 




min; 


ni.mi 


ni.mi 


r = a + km 


or 



where (a Xy a y ,a z y is any point on the line and (m Xl m yy m z y is the direction of the line. Without any loss of 
generality, we could set a z = 0.0 and m z = 1.0 (al though this is not true for lines that are parallel to the z- 
direction, we assume that the through plane direction is properly oriented such that no line is parallel to the 
z-plane). Thus the equation for a line is 





For a particular z, we estimate the (x, y) component of the line. By repeating the procedure for a series of 
(zi,...,z„), we obtain (n,... ,x„) and (y u . ..,}/„). We could form the following equation for the x components 


f *i ') 

\ Xn / 


/ 1 Zi 


\ 1 Zn 



Cz, 1 


(14) 


where (e x ,i , . . . , e x>n y are the estimation errors for the (xi , . . . , x n )'. Based on the assumption that the estimation 
errors for each z are independent and have the same variance (which is true in most cases), (a Zy m x ) in (14) could 
be estimated by any standard linear regression technique [23]. Using the following notation: 


x — 





the estimated values for ( a xy m x ) are 


and 



) = (Z T Z) _1 Z T x 


=var(a x ) = cr 2 x {( Z^Z)- 1 ^ 
a m,x = var(m x ) = a\{{ Z T Z ) -1 } 2 2 


(15) 

(16) 


where a\ = uar(e<) and {(Z T Z)- 1 } y refers to the (i,j) element of the matrbc (Z T Z) _1 . We perform the same 
procedure for the y component of the line. 


F Estimation of a 2 

In cases where var(e) is unknown, we could use the following method to obtain an estimated a 2 [23]. Since a 2 
is the expected squared value of an error, e*, it is natural to use the simple squared value of the residuals. The 
vector of residual for the x component is 


x — X 

x — Z(Z T Z) _1 Z T x 


Under the assumption that the errors axe uncorrelated with constant variance a 2 , an unbiased estimate of a 2 is 


* l = 


n-p 


where p is the number of parameter (which is 2 in this case) and n is the data size, cr 2 is estimated in a similar 


manner. 
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Registration of Serial Skeletal Images for Accurately Measuring 

Changes in Bone Density 

Chye Hwang Yan, G.S. Beaupre *, R.T. Whalen**, S. Napel *** 

Department of Electrical Engineering, Stanford University, CA 94305, USA 


Relevance to Musculoskeletal Condition: The non-invasive determination of changes in bone mass is 

valuable as a clinical assessment tool in a variety of settings. Examples include: evaluating the efficacy of drug 
therapy in the treatment of osteoporosis; understanding local bone adaptation in response to exercise intervention; 
and assessing the progression of bone loss with skeletal disease. 

Introduction: The accurate measurement of bone apparent density from Computed Tomography (CT) requires 
a consideration of image noise, marrow presence and beam hardening. Non-invasive measurement of changes in 
bone apparent density (grams/cm3) further necessitates imaging the same skeletal site at different times and 
registering the serial images. Mis-registration can introduce errors, particularly when density gradients are high. 
The effect of mis-registration errors generally increases as the size of the volume of interest (VOI) decreases. For 
this reason and because of image noise relatively large VOI's (3000-4000 mm3) are typically analyzed [1,2]. 

We have developed a semi-automated, 3D, surface-based registration technique that does not require the use ofan 
external frame or fiducial markers. The accuracy of our non-invasive registration technique is comparable to that 
which can be obtained using standard invasive approaches (e.g., skeletal pins and frames). The objective of this 
study was to assess the accuracy and precision of repeated registration of the human calcaneus whe n imaged 
with a clinical CT scanner. 

Methods: Excised human calcanei were secured within a custom-made fiducial frame, referred to as a ” tiger 

cage” (Fig. 1). The tiger cage/calcaneus was scanned in different orientations using a GE HiSpeed Advantage 
CT scanner. Scanning technique consisted of: helical mode; 120 kVp; 240 mA; 1.0 mm slice thickness; pitch = 
1.7, scan field of view = 250 mm; display field of view = 120P130 mm; standard reconstruction algorithm. Using 
a 512 x 512 reconstruction matrix, the in-plane pixel size was approximately 0.25 mm. Reconstructions were 
done every 0.5 mm in the through-plane direction. To examine the effects of image noise and contrast, the tiger 
cage/calcaneus was scanned both in air and in water contained within a 240 mm diameter, thin-walled (12.7 mm) 
plexiglass cylinder. 

We have recently developed an improved gold standard [3] that we used to assess the accuracy of our new surface- 
based registration of the calcaneus. We determined registration accuracy by: 1) obtaining a gold standard 
transformation for the tiger cage between pairs of scans; 2) obtaining a transformation for the calcaneus using 
surface points only (200,000 to 300,000 points per calcaneus); 3) transforming all calcaneus voxels using both 
the gold standard and the calcaneus transformations; 4) calculating the distance (error) between corresponding 
voxels in the two transformed calcaneus images. 

Results: Shown in Table 1 are the registration errors for the calcaneus scanned in air and scanned in water. 
Note that the average errors in air and water are similar and sub-pixel in magnitude. 

Discussion: We have developed a highly accurate frameless registration technique for use in determining 

changes in bone density. The average error throughout the volume is sub-pixel in magnitude. It should be 
noted, however, that our registration errors are spatially correlated; highest accuracy exists near the centrum of 
the calcaneus. To our knowledge this new technique permits accurate registration of bone volumes an order of 
magnitude smaller than previously possible using clinical CT scanners. The technique is robust in the presence 



Figure 1: Calcaneus within tiger cage. 
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of noise with minimal degradation of registration accuracy when scanning with a water bolus. Furthermore, by 
using surface points only, the accuracy of registration is independent of changes in bone apparent density. 

Highly accurate and precise registration is an essential component for non-invasive monitoring of changes in local 
volumetric bone density. Using this newly developed registration technique it will be possible, for example, to 
monitor more accurately osteoporosis progression and treatment in both large and small regions using a clinical 
CT scanner. In particular, the ability to measure bone changes in small regions allows one to examine relationships 
between remodeling rates and local bone morphological parameters such as apparent density and related variables 
such as bone specific surface area. 

References: [1] Klotz et al. Trans. Med. Imaging, 37IP376, 1989; [2] Tanno et al. Bone, 239P247, 1996; [3] 
Yan et al., accepted 1996 Radiological Soc. N. Amer. 
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Precise and Accurate Gold Standard for Multimodality and Serial 

Registration Method Evaluations 
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Purpose: The accuracy of a registration method is usually determined by using the registration of an attached 

stereotactic frame as the gold standard. Existing frame registration techniques are: (1) matching each 2D slice 
separately to an N-bar model (2) using 2 detected N-bars (6 per frame) as the coordinate system. These techniques 
assume a precise frame model and only utilize partial 3D information. We propose a technique that assumes only 
linear structures (rods) in the frame and utilizes all 3D information. 

Materials and Methods: Our technique fits lines to all linear structures in each data set and registers them 

using a closed form solution. We obtained CT scans (0.24x0. 24x1. 00mm 2 3 voxels) of a frame at 8 orientations. We 
performed 28 registrations between scans using existing and proposed techniques. We defined registration error 
as the maximum deviation between registered corresponding rods. 

Results: Registration errors (average, maximum) for the existing 2D and 3D techniques were (0.18mm, 0.22mm) 
and (0.08mm, 0.11mm) respectively. The proposed method reduced the error to (0.035mm, 0.04mm). 

Conclusion: The proposed technique produces more accurate registration without requiring any knowledge of 

the frame geometry. It is used to evaluate our new surface registration technique, for high-resolution QCT, which 
has an accuracy comparable to existing gold standards. 


Take Home Points: (1) accurate registration with simple frame. (2) global minimum registration. 
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